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EXECUTIVE  SUMMARY 


In  this  work  we  further  develop  a  linear  time-varying  system  model  of  sea  clutter  returns  begun  in  a 
previous  NRL  report  [1],  We  present  analytic  formulae  for  the  clutter  autocorrelation  function,  which  can 
be  used  in  the  design  and  analysis  of  clutter  mitigation  signal  processing  techniques.  We  further  refine  the 
model  to  account  for  the  periodic  nature  of  the  pulse  ambiguity  function  and  refine  the  computation  of  the 
Doppler  spectrum  to  account  for  wave  motion  as  well  as  improve  its  numerical  integrability. 

We  apply  these  tools  to  demonstrate  how  the  effectiveness  of  a  candidate  waveform  in  the  presence  of 
clutter  may  be  evaluated. 


E-l 


APPLICATION  OF  A  STATISTICAL  LINEAR  TIME-VARYING  SYSTEM  MODEL  OF 
HIGH  GRAZING  ANGLE  SEA  CLUTTER  FOR  COMPUTING  INTERFERENCE 

POWER 


1.  INTRODUCTION 

Statistical  linear  time-varying  (LTV)  system  theory,  particularly  in  the  case  of  wide-sense  stationary 
uncorrelated  scattering  (WSSUS)  was  first  put  on  firm  mathematical  footing  by  Bello  [2],  but  was  hinted 
at  by  other  authors  that  preceded  him,  such  as  Price  and  Green  [3,  4],  WSSUS  channel  models  have  been 
widely  used  in  communications  [5-7],  but  have  only  occasionally  seen  use  in  radar  circles  [8-10]  despite 
originating  largely  in  this  community.  In  this  work  we  further  refine  a  WSSUS  model  for  a  sea  clutter 
scattering  channel  developed  by  the  authors  [1]  and  demonstrate  how  it  can  be  used  to  evaluate  waveform 
performance  in  the  presence  of  clutter. 

Let  the  single -polarization  response  y(t)  of  a  monostatic  radar  scattering  channel  to  input  signal  x(t )  be 
given  by  the  superposition  integral: 


/OO 

h(r,t)x(t  —  r)dr  (1) 

-OO 

where  h(r,  t )  is  the  time- varying  impulse  response.  Note  that  if  there  is  no  relative  motion  in  the  channel 
between  the  radar  and  the  scatterers,  h(r.  t )  =  h(r),  and  the  system  becomes  a  linear  time-invariant  (LTI) 
system  whose  output  is  given  by  the  convolution  of  x  and  h. 

In  applying  this  model  to  a  radar  system,  we  use  the  following  assumptions: 


•  The  LTV  impulse  response  h(r,  t)  is  only  valid  over  a  single  coherent  processing  interval  (CPI) 

•  Over  a  single  CPI,  the  range  to  target(s)  is  roughly  constant 

•  The  lag  variable  r  represents  downrange  delay 

•  The  time  variable  t  represents  the  change  in  the  channel  over  time 


Variation  in  h  with  respect  to  t  represents  fluctuation  in  the  signal  amplitude  due  to  relative  motion  between 
the  radar  and  the  scatterers  and,  subsequently,  produces  a  Doppler  shift  in  the  received  signal  that  is  not 
present  if  the  system  is  LTI. 

The  other  key  assumption  of  a  WSSUS  model  is  that  the  channel  is  random,  and  it  is  usually  assumed 
to  have  Gaussian  statistics  in  both  the  in-phase  and  quadrature  channels.  If  the  channel  consists  only  of 
clutter,  this  corresponds  to  the  Rayleigh  amplitude  model  commonly  used  for  high-grazing  angle  clutter. 
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which  result  derives  from  the  assumption  that  a  large  number  of  independent  scatterers  are  contained  in  a 
single  resolution  cell  [11]. 

Note  that  due  to  the  lineality  of  (1),  if  multiple  scatterers  are  present,  the  overall  impulse  response  is 
given  by: 


h{r,  t)  =  hi(r,  t )  +  h2{r,  t) + 

where  h,  represent  the  impulse  responses  of  individual  scatterers. 


(2) 


2,  AUTOCORRELATION  FUNCTION 

Knowledge  of  the  clutter  autocorrelation  function  allows  one  to  predict  the  effectiveness  of  signal  pro¬ 
cessing  techniques  such  as  space-time  array  processing  (STAP)  [12],  therefore  being  able  to  model  it  effec¬ 
tively  will  be  invaluable  in  the  analysis  and  design  of  future  signal  processing  algorithms.  The  autocorrela¬ 
tion  function  of  the  random  channel  h(r,  t)  is  given  by  [5]: 

E  \h*(r,  t)h{r' ,  t  +  At)]  =  A(t,  At)S(r  —  t'),  (3) 

which  implicitly  states  that  the  channel  is  wide-sense  stationary  (WSS)  in  t,  and  it  undergoes  uncorrelated 
scattering  (US)  in  r,  hence  the  name  WSSUS.  The  function  A(t,  At)  is  known  simply  as  the  WSSUS 
autocorrelation  function. 

The  scattering  function,  S(r,  p),  which  is  a  representation  of  the  channel  gain  as  a  function  of  delay  and 
Doppler  shift,  is  simply  the  Fourier  transform  of  A(t,  At)  along  the  At  axis: 


/OO 

A(t,  At)e-j2npAtdAt. 

-OO 


(4) 


Note  that  if  multiple  scatterers  are  present,  as  in  the  case  of  (2),  and  if  those  scatterers  are  statistically 
independent,  it  can  be  shown  that  the  overall  scattering  function  is  given  by: 


S(r,p)  =  Si{r,p)  +  S2{t,p)  +  . . . , 


(5) 


where  S)  represent  the  scattering  functions  of  the  individual  scatterers.  This  superposition  property  also 
applies  to  the  autocorrelation  functions  (3),  provided  the  target  responses  are  independent. 

In  the  author’s  previous  work  [1],  it  was  shown  that  the  WSSUS  autocorrelation  function  A(t,  At)  for 
a  sea  clutter  scattering  channel  with  no  other  targets  is  given  by: 


\2(T°C 

(47t)3(4ct)3 


G2(<M)ei2^pW)Ai#, 


A(t,  At) 


(6) 
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where 


p(<p,  6)  =  px  cos  9  cos  <p  +  pz  sin#  (7) 

is  the  observed  Doppler  shift  as  a  function  of  azimuth  angle  <l>  and  depression  angle  9  for  an  airborne  radar 
moving  at  velocity  v  =  [vx,  0,  vz]T,  A  is  the  carrier  wavelength,  c  is  the  speed  of  light,  px  =  2vx/\, 
pz  =  2vz/\,  and  a0  =  a0(a(r ))  is  the  normalized  radar  cross  section  (NRCS)  of  the  clutter  as  a  function 
of  grazing  angle  a,  which  itself  is  a  function  of  r.  We  will  assume  an  antenna  power  pattern  G{d>.  6)  given 
by: 


G(4>,  9)  =  D 


sine 


sine 


0 e  -  9a ) 


(8) 


where  D  is  the  antenna  directivity,  W  =  1.20671  rad  is  a  normalization  constant,  B  is  the  3  dB  beamwidth 
in  radians,  and  {(bn.  9a)  are  the  orientation  angles  of  the  main  beam.  We  can  approximate  one  of  the  sine 
factors  using  the  Dirichlet  kernel  to  facilitate  computation  of  the  integral  in  (6)  as  follows: 


sine 


O  -  <l>a 


sin(7T W(<p  -  (pa) / B) 

( 2nW/B )  sin ((</>  -  <pa)/ 2) 


(9) 


If  we  define  the  quantity  N: 


W 

n  =  2*b 


(10) 


then  we  can  express  (9)  as  follows: 


W 

sine  I  —  (0  -  (pa 


sin (N(cp  -  (pa)/ 2) 


N  sin((</>  -  (pa)/ 2) 


(11) 


which,  if  N  3>  1,  then  to  good  approximation  we  can  round  N  to  the  nearest  integer  and  express  (1 1)  by: 


sin (N((p  -  (pa)/2) 
N  sin((0  -  (pa)/ 2) 


1 

N 


N- 1 


e3m,{<j>-<j>a) 

m= 0 


(12) 


Now  we  can  solve  the  integral  in  (6)  as  follows: 


G2((p,9)ej27rp'^'9)Atdd) 


W 

~B 


(' 0  ~  9a) 


^ j2npz  sin  6 At 


2 

ej2irpx  cos  6  cos  <j>A 


=  D2  sine2 


(13) 
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Note  that  if  we  expand  the  first  factor  in  the  integrand  in  a  Fourier  series: 


TV— 1  z 

_L  pjm(<t>-<t>a) 

ivZ. 


y  c^jk^ 


where 


(F-§)e_^“>  \k\<N 


and  then  substitute  the  Jacobi- Anger  expansion  [13]  for  the  second  factor  of  the  integrand  in  (13): 


0j(2'Kpx  cos  6  At)  cos  (ft  _ 


y,  jnJn{2irpx  cos  9At)ejn<l>, 


where  Jn{x)  is  the  nlh -order  Bessel  function  of  the  first  kind,  we  can  solve  the  integral  (13)  as  follows: 

N—l  2 

[  1  y  epn(<p-4>a)  ej2npxCoS0coS<t>Atd(j>=  f  y  CfcejW  y  -n  jn^Px  cos  9 At)ejn<t,d(j) 

J{2*)  N  ^  n 

=  y  y  ckjnjn(2irpx  cos  9 At)  [  ej(k+n)4>d<j) 

k  n  W 

=  EE  CkjnJn(2irpx  cos  9At)2ir5[k  +  n] 

k  n 

=  27T  y  C—njn  J n  (27T px  COS  9  At) 
n 

N  1  1  \ 

=  2vr  E  (  N  ~  W  )  e3n4>a3njn{^px  cos  9 At) 


N  /  I  IN 

E  (  1  -  -J7  )  ej^a+^nJn(27Tpx  cos  9 At) 

_  Ar  \  / 


The  correlation  function  (6)  is  then  given  by: 


\2(t°cttD2 

(4tt)3(^ct)3N 


^(0-  9a)^j  ■  ei2'Kp~ sin eAt  ■  y  (i  -  M)  ej^a+n^nJn{2Trpxcos9At). 


(18) 
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Fig.  1 :  Unsmoothed  Clarke  spectrum  with  G  =  1,  px  cos  9  =  1  Hz,  and  pz  sin  6  =  0  Hz. 


3,  DOPPLER  SMOOTHING 


3.1  Frequency  Domain 


From  the  author’s  previous  report  [1],  it  was  determined  that  the  clutter  scattering  function  S(t,  p)  was 
given  by: 


S(t,  p) 


(  A 2c  .  ar°(a(r))  _  _ G,2(</>(r,p),6>(r)) _ 

(4FF  (|cr)3  ^/(Px  COS(0(T)))2 -(p-pz  sin (0(r)))2  ’ 


sin(6»(r))|  <  px  cos (0(r)) 

•  (19) 
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else 


This  is  a  generalization  of  the  famous  result  given  by  Clarke  in  his  landmark  paper  [14].  However,  this 
result  was  derived  assuming  the  antenna  pattern  Gif.  9)  was  an  even  function  of  ([>.  For  the  most  general 
case,  i.e.  when  cj)a  /  0,  the  scattering  function  is  given  by: 


S(t,  p) 


f  A 2Cff°(a(r))  _  1  G2{(P(t,p),6(t))+G2(-<P(t,p),0(t)) 
(47t)3( |cr)3  2  y/ (px  cos {9{t)))2 -{p- pz  sin (6(t)))2  ’ 


\p  -  pzsin(0(r))|  <  pxcos{9(t)) 


(20) 


0 


else 
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Fig.  2:  Doppler  spectrum  smoothing  kernel. 


In  this  section,  we  will  define  the  function  F(p)  to  represent  the  factors  of  (20)  that  only  depend  on  p: 

F{P) 

This  function  is  plotted  in  Figure  1.  It  should  be  noted  that  (21)  assumes  the  clutter  is  not  moving.  As  well, 
it  should  be  noted  that  (21)  possesses  singularities  when  \p  —  pz  sin  |  =  px  cos  9,  which  presents  problems 
during  numerical  integration  because  it  makes  the  result  very  sensitive  to  the  sample  spacing. 

One  way  of  alleviating  both  of  these  problems  is  to  smooth  the  Doppler  spectrum  by  convolving  the 
scattering  function  (20)  with  a  smoothing  kernel  I\(p)  to  remove  the  singularities  and  broaden  the  Doppler 
spectrum  by  a  modest  amount  to  attempt  to  model  wave  motion.  In  this  work,  we  choose  a  parabolic 
smoothing  kernel  given  by: 


g2(Mp),0) 


sj p*.  cos2  9  —  (p  —  pz  sin  8)2 


(21) 


K(p) 


4  pm  \  2  pm  J 


(22) 


This  function  was  chosen  because  it  has  a  mostly  “round”  shape  and  (more  importantly)  the  convolution  of 
K  with  F  possesses  an  analytic  solution.  Note  that  K(p)  is  normalized  so  that  it  integrates  to  one  so  that 
the  total  power  is  conserved.  The  Doppler-domain  convolution  (of  the  factors  that  depend  on  p)  is  expressed 
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This  result  is  plotted  in  Figure  3.  Note  that  now  the  singularities  have  been  removed  and  the  Doppler 
spectrum  has  been  broadeded,  as  desired. 
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Fig.  3:  Smoothed  Clarke  spectrum  with  G  =  1,  pm  =  0.15  Hz,  px  cos  9  =  1  Hz,  and  pz  sin  8  =  0  Hz. 
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Fig.  4:  Time-domain  autocorrelation  window  function  used  to  account  for  wave  motion. 


3.2  Time-Domain 

To  obtain  an  expression  for  the  time-domain  windowed  autocorrelation  function  Aw(r,At),  we  will 
perform  Fourier  inversion  on  the  smoothed  clutter  scattering  function  Sc(r,  p)  to  yield  the  time-domain 
autocorrelation: 


Sc(t,  p)  =  S (t,  p)  *K(p )  (27) 

=>  A w(t,  At)  =  F~x  [Sc(r,  p)\ 

=  f~x  [S(T,  p)  *  K(p)\  =  F-x  [S(T,  p)]  f;1  [k(p)} 

=  A(t,  At)k(At),  (28) 


where 


k(At)  =  F-1  [. K(p )]  =  3 


/  sin(27rpmAi) 
V  (27T/3mAf)3 


COs(27T Pm,  At)  \ 

(2irpmAt)2  ) 


(29) 


Note  that  pm  can  be  modified  to  include  range-dependence,  i.e.  prn  =  pm{r)  without  loss  of  generality. 
The  window  function  k(At)  is  plotted  in  Figure  4.  The  resultant  autocorrelation  can  then  be  found  by 
substituting  (18)  into  (28). 

The  Python  code  used  to  generate  Figures  1-4  is  found  in  Appendix  A. 
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4.  OUTPUT  TIME-FREQUENCY  POWER  DISTRIBUTION 

We  will  assume  the  matched-filter  output  is  the  correlation  of  N  pulses  against  an  infinite  pulse  train 
input,  represented  by  the  periodic  ambiguity  function  (PAF)  \xnt\-  given  by  [15]: 


I Xnt{t,  p)| 


1 

NTr 


u(t)u*(t  +  T)ei27Tptdt 


(30) 


where  u(t)  is  the  normalized  unit-energy  input  pulse  train  and  Tr  is  the  pulse  repetition  interval  (PRI).  The 
pulse  repetition  frequency  (PRF)  is  given  by  Fr  =  1/Tr.  The  PAF  has  the  following  properties  that  arc 
relevant  to  our  application: 


\xnt(t,p)\  =  \xnt(t  +  Tr,p)\ 


\xnt(t,  p) | 


p)  | 


sin  N7rpTr 
N  sin  7r  pTr 


The  time-frequency  power  distribution  at  the  channel  output  is  given  by: 

P(r,p)  =  ExS(r,p )  *  *\xnt{t,  p)\2 . 

=  EXJ  J  S(t',p')\xnt(t  ~t',p-  p')\2dT'dp'i 


(31) 

(32) 


(33) 


where  Ex  is  the  energy  of  the  pulse  train.  Note  that  since  \xnt\2  is  periodic  along  the  r  axis,  the  result  of 
the  convolution  integral  (33)  is  periodic  along  this  axis  as  well.  If  we  assume  the  support  of  S  is  limited  to 
[0,  LTr\  for  some  positive  integer  L  along  the  r  axis,  then  the  output  power  distribution  can  be  written  as  a 
circular  convolution: 


P(t,  p)  =  Ex 


pLTr  coo 

/  /  S(t',p')\xnt{t  -  r',p  -  p')\2dp’dT' 

J  0  J  —  oo 

=  ExS(t,  p)®*  \xnt{t,  p) I2, 

T  P 


(34) 


where  ©  represents  circular  convolution  along  the  r  axis  and  *  represents  ordinary  convolution  along  the 

r  p 

p  axis.  The  problem  was  cast  into  this  form  so  that  discrete  implementations  can  use  the  Fast  Fourier 
Transform  (FFT)  to  efficiently  compute  the  periodic  output  power  distribution  along  the  r  axis  without 
zero-padding. 
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5.  SIGNAL-TO-CLUTTER  RATIO 

We  will  now  use  the  results  from  the  previous  sections  to  compute  target  signal  and  clutter  signal  power 
distributions.  We  will  use  as  our  transmitted  signal  an  unmodulated  train  of  N  =  8  rectangular  pulses  with 
pulse  width  T  =  10  /as,  PRI  Tr  =  100  /is,  and  Ex  =  N  J.  A  plot  of  the  PAF  for  this  pulse  train  is  shown 
in  Figure  5  for  the  principal  PRI  and  PRF  interval.  A  sample  rate  Ts  =  T / 10  =  1  /rs  and  a  frequency  bin 
spacing  of  A p  =  Fr/ 256  =  39.0625  Hz  was  used  for  all  plots. 

A  smoothed  clutter  scattering  function  Sc(t.  p )  was  generated  from  (20)  and  (27)  with  the  following 
parameters: 


/  = 

10 

GHz 

V  = 

[vx 

,Vy,Vz]T 

=  [300,  i 

Pm  1 

=  66.713  Hz 

G(4> 

,0) 

given  by 

(8)  with: 

- 

D 

=  100 

- 

B 

=  10° 

- 

<t>a 

=  0° 

— 

0a 

=  25° 

Fig.  5:  Periodic  ambiguity  function  (PAF)  of  a  train  of  N  =  8  rectangular  pulses. 
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Pc(t,  p )  (dB  normalized  to  peak) 


0  ^ 
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Fig.  6:  Time-frequency  distribution  of  clutter  signal  power. 


•  cr°(a)  data  taken  from  Figure  7.13  of  [16]. 


with  the  following  earth  geometry  assumptions: 


•  h  =  5  km 

•  a(r)  =  0(t)  =  sin_1(/i/(^cT))  ==>  flat  earth  model. 

•  p )  =  cos"1  ),  assuming  0  <  (j)  <  vr. 


The  support  of  ,Sc(r:  p)  along  the  r  axis  was  assumed  to  be  [0,  500  x  10~6]  s,  implying  that  L  =  5.  The 
relation  derived  in  (34)  was  used  to  compute  the  clutter  distribution,  and  the  result  Pc(r,  p)  is  plotted  in 
Figure  6.  From  this  figure  the  altitude  return  at  r  «  33.4  //s  is  clearly  visible  as  well  as  the  mainlobe  clutter 
at  t  «  78.9  /is  and  p  «  8.14  kHz. 

The  time -frequency  distribution  of  a  point  target  with  RCS  a  =  40  dBsm  located  at  (r,  p)  is  given  by 

[1]: 


or  \  o  G'2(^>('r’  6'(t))A2ct  _  2 

Pt{T,p)  =  Ex - j - |xiVT(r-r,p-p)|  , 


(47T 


(35) 
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0  20  40  60  80 

T(flS) 


Fig.  7 :  Time-frequency  distribution  of  target  signal  power. 


which  is  simply  a  shifted  and  scaled  version  of  the  PAF.  If  the  main  beam  is  pointed  at  the  target  then 
f  ~  78.93  fi s  and  p  =  18.14  kHz.  However,  due  to  the  fact  that  the  digital  receiver  can  only  observe 
spectrum  from  0  to  Fr,  the  main  lobe  of  the  power  distribution  will  appeal-  at  p  —  Fr  ~  8.14  kHz.  The 
result  is  plotted  in  Figure  7.  The  signal-to-clutter  ratio  (SCR)  Pt/Pc  as  a  function  of  (r,  p)  is  calculated  and 
plotted  in  Figure  8.  The  observed  SCR  in  the  direction  of  the  main  lobe  was  computed  as  a  function  of  N 
and  is  shown  in  Table  1 . 


N 

SCR  (dB) 

1 

7.228 

2 

7.331 

4 

7.463 

8 

7.730 

16 

8.414 

Table  1 :  Signal-to-clutter  ratio  as  a  function  of  num¬ 
ber  of  pulses  coherently  integrated. 
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Fig.  8:  Signal-to-clutter  ratio  as  a  function  of  delay  and  Doppler  shift. 


It  can  be  seen  from  Table  1  that  increasing  the  number  of  pulses  integrated  does  improve  the  SCR,  albeit 
modestly.  If  the  clutter  return  from  each  pulse  was  completely  independent,  the  SCR  would  be  expected  to 
increase  linearly  with  N.  However,  because  the  pulses  arc  highly  correlated,  integrating  extra  pulses  does 
not  improve  the  SCR  very  much. 

The  Python  code  used  to  generate  Figures  5-8  is  found  in  Appendix  B. 
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6.  CONCLUSION  AND  FUTURE  WORK 

In  this  report,  further  derivations  related  to  the  implementation  of  a  WSSUS  clutter  model  have  been 
furnished,  and  output  time-frequency  power  distributions  have  been  produced  to  simulate  what  an  actual 
radar  would  “see”.  Python  code  for  said  simulations  has  also  been  provided. 

Future  work  on  this  topic  could  focus  on  verifying  the  assumptions  behind  the  WSSUS  model  against 
experimental  data  and  applying  this  model  to  the  analysis  and  design  of  practical  clutter  mitigation  tech¬ 
niques. 
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Appendix  A 

CODE  LISTING:  FIGS1.PY 


from  pylab  import  * 

close ( ’  all  ’  ) 

from  figs2  import  doppler_slice 

def  main  ( ) : 

N  =  2**14 

dt  =  lin sp ac e  (  —  8 . ,  8.,  N) 
rho  m  =  1. 

figsize  =  (5.25,3.5) 


def  k ( 1 1  ) : 

x  =  2*  pi  *rho_m*  1 1 

rval  =  3.*(sin(x)/x**3.  —  cos(x)/x**2.) 

r  v  a  1  [  x  ==  0 .  ]  =  1 . 
return  rval 

figure  (  figsize  =  figsize  ) 
plot ( dt ,  k(dt),  ’k’) 
xlabel  (r  ’  $  \  rho  _m  \  Delta^t$  '  ) 
ylabel(r'$k(\Delta^t)$’) 
ylim  (  — 0.2  ,  1.1) 

grid  () 

t  ig  ht  _1  ay  o  u  t  () 

savefig  (  ’  fig_kdelta  .pdf’) 


def  K(ff): 

rval  =  3  . /( 4  .  *  rho_m ) * ( 1 .  —  (  f f  / rho_m  ) *  * 2  . ) 
rval[abs(ff)  >=  rho  m]  =  0. 
return  rval 

rho  =  li  nsp  ac  e  (  —  1 .5  ,  1.5,  N) 

figure  (  figsize  =  figsize  ) 
plot  (rho,  rho_m*K(  rho  )  ,  ’k’) 

xlabel ( r ’ $\rho /\ rho_m$  ’ ) 
y  label  (r  '  $\rho_m  JC(\ rho  )  $  ’  ) 
xlim  (min( rho  )  ,  max(rho)) 
ylim  (—0.1,  0.8) 
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grid  () 

tight_layout  () 

savefig  (  ’  fig.krho  .  pdf’  ) 


rho_x  =  1. 
rho  m  =  0.15 

theta  =  zero s ( rho . shape  ) 
rho_z  =  0. 

def  unsmoothed_slice  (  rr  ) : 

i  =  abs  (rho— rho_z  *  sin  ( theta  ))  <  rho_x*cos  ( theta  ) 
returnval  =  zero s ( rr . shape ) 

returnval[i]  =  l./sqrt((rho_x*cos(theta[i]))**2.  —  \ 
(rr  [i]  —  rho_z*sin(theta  [i  ] ) )  *  *  2  . ) 
return  returnval 

figure  (  figsize  =  figsize  ) 

plot(rho,  unsmoothed_slice  (rho  )  ,  ’k’) 

xlim  (min( rho  )  ,  max(rho)) 

ylim  (  — 0.5  ,  4) 

xlabel (r ’$\rho$^(Hz) ’ ) 

ylabel (r ’ $F(\ rho ) $  ’ ) 

grid  () 

tight.layout  () 

savefig  (  ’  fig_f_unsmoothed  .  pdf  ’ ) 

F  =  doppler .slice  (rho  ,  rho_x  ,  rho_z  ,  rho_m  ,  theta  =  theta) 

figure  (  figsize  =  figsize  ) 

plot (rho  ,  F,  ’k’ ) 

xlim  (min( rho  )  ,  max(rho)) 

ylim  (  — 0.5  ,  4) 

xlabel ( r ’ $\rho$  ^(Hz) ’ ) 

ylabel ( r ’ $\{F*K\  }(\rho)$’) 

grid  () 

tight.layout  () 

savefig  (  ’  fig_f_smoothed  .pdf’) 
show  () 


if  name  ==  ’  main  ’  : 
main  () 


Appendix  B 

CODE  LISTING:  FIGS2.PY 

from  pylab  import  * 

close (  ’  all  ’  ) 

SPEED  _OF_LIGHT_M_PER_S  =  29979245  8.0 

#  Helper  functions 
def  db(x): 

return  10*logl0(x  +  le— 99) 

def  cosd(x): 

return  cos ( x* pi / 1 8  0  . ) 

def  sind (x  ) : 

return  sin  ( x* pi / 1 8 0 . ) 

def  tand(x): 

return  sind (x  )/ cosd (x) 

def  esc (x ) : 

return  1  ./  sin  (x) 

def  cscd (x  ) : 

return  c sc ( x* pi / 1 8  0  . ) 

def  nnz (x ) : 

#  returns  number  of  nonzero  elements 
return  sum  ( ( x  !  =  0 ) .  f  1  a  1 1  e  n  ( ) ) 

def  circconv2_tauaxis  (x ,  y): 

#  circular  convolution  along  second  dimension 

#  normal  zero— padded  linear  convolution  along  first  dimension 
nix ,  n2x  =  x . shape 

nly ,  n2y  =  y . shape 

if  n2x  !=  n2y  :  raise  Exception  (’ second  ^dimensions  ^must^be  ^equal  ’ ) 
nl  =  nix  +  nly  —  1 
npadx  =  nl  —  nix 
npady  =  nl  —  nly 
xpad  =  concatenate  (( 
complexl28  (x)  , 
zero s (( npadx ,  n2x), 
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dtype=complex  , 

))  ,  axis  =0) 
ypad  =  concatenate  (( 
complexl28  (y )  , 
zero s ( ( npady ,  n2y), 
dtype=complex  , 

))  ,  axis  =0) 

return  i  f f  1 2  (fft2(xpad)*fft2(ypad)) 

def  si  nc  2p  at  t  ern  (  phi.deg  ,  theta.deg  ,  gain_db=20,  beamwidth.deg  =  1 0) : 

#  sine  "2  power  pattern  in  radial  direction 
phi  =  phi  _deg  *  pi  /  1  8  0 . 

theta  =  theta_deg*pi/180. 
gain  =  1 0 . *  * ( g ain_db / 1  0  . ) 
beamwidth  =  beamwidth_deg*pi /I  80. 

W  =  0.44294647068945234030836999*2  #  sine  half —power  width 

return  gain*  sine  (  sqrt  (  phi  **2.  +  t he t a  *  *2 . ) *W/ beamwidth ) * *2 

def  sinepattern  (  phi_deg  ,  theta.deg  ,  gain_db=20,  beamwidth.deg  =  1 0) : 

#  sine  power  pattern  in  each  direction 
phi  =  phi  _deg  *  pi  /  1  8  0 . 

theta  =  theta_deg*pi/180. 
gain  =  1 0 .  *  *  (  g  ain_db  /  1  0  . ) 
beamwidth  =  beamwidth.deg  *  pi  /  1  8  0 . 

W  =  0.6033545  64401614197645753  84*2  #  sine  half —power  width 

return  gain  *abs  (  sine  (  phi  *W/  beamwidth  )*sinc(theta  *W/  beamwidth  ) ) 


def  plot_image  (x  ,  y  ,  z  ,  f  i g  s i  z e  =None  ) : 
dxl  =  x[l,0]  —  x[0,0] 
dx2  =  x[0,l]  —  x[0,0] 
dx  =  max(dxl,dx2) 

dyl  =  y  [  1  ,0]  -  y  [0  ,0] 
dy2  =  y  [0 , 1  ]  -  y[0,0] 
dy  =  max(dyl,dy2) 

if  figsize  !=  None: 

figure  (  figsize  =  figsize  ) 
else  : 

figure  () 
imshow  ( 
z  , 

aspect= ’ auto ’  , 
origin=’lower’  , 
i  n t  e  rp  o  1  a t  i  o  n  =  ’  none  ’  , 
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extent  =( 

x .  min  ()  —  dx /2 .  , 
x  .  max  ()  +  dx  /  2  .  , 

y .  min  ()  —  dy /2 .  , 
y  .  max  ( )  +  dy  /  2  .  , 

), 

cmap=cm.  gray  , 

) 

colorbar  () 

def  p  ul  s  e  _  tr  ai  n  _f  ac  tor  ( rho  ,  N,  p,  T_r): 
return.val  =  zeros ( rho . shape ) 
i  =  (rho  ==  0) 

r  etur  n  val  [~  i  ]  =  abs  (  sin  (  pi  *  rho  [~  i  ]  *  (N—  abs  (p ))  *  T  r  )  /  sin  ( 
pi  *rho  [~  i  ]*  T.r  )) 
return.val [i]  =  N  —  abs (p) 
return  return.val 

def  doppler  .integral  (rho  ,  L,  rhox  ,  rhoz  ,  rhom  ,  theta): 
rhox  +=  le— 8 

sqrtarg  =  rhox**2.*cos(theta)**2.  —  (rhoz*sin(t  he  t  a )— L)  *  *  2 . 

badinds  =  sqrtarg  <  0 

sqrtarg [ badinds ]  =  0. 

sqrtpart  =  sqrt(sqrtarg) 

arctanpartargnum  =  rhoz*  sin  ( theta)— L 

arctanpart  =  zeros ( badinds . shape ) 

arctanpart  [badinds]  =  sign  (arctanpar  tar  gnu  m  [  badinds])*  pi  /  2. 
arctanpart  [~  badinds  ]  =  arctan  ( 

arctanpartargnum  [~  badinds]  /  sqrtpart  [“badinds]) 
return  3  .  /(  8  .  *  rhom  *  *  3  . )  *  ( 

( rhox * *2. *  cos ( t he t a ) * *2 .  +  2 . * rhoz * *2 . *  sin ( t he t a ) *  *  2 .  —  \ 
4*rhoz*rho*  sin  ( theta  )  +  2.*rho**2.  —  2.*rhom**2.  \ 

)*  arctanpart  +  s  qr  tp  ar  t  *  (3 .  *  rhoz  *  sin  ( the  t  a  )  — 4.  *  rho  +  L)) 

def  doppler  .slice  (rho  ,  rhox,  rhoz,  rhom,  theta): 

return.val  =  zeros ( rho . shape ) 

11  =  abs(rho  —  rhoz* sin ( theta  ))  <=  rhox*cos ( theta )  —  rhom 

12  =  abs(rho  —  rhoz* sin ( theta )  +  rhox *  cos ( t he t a ) )  <  rhom 

13  =  abs(rho  —  rhoz* sin ( theta )  —  rhox *  cos ( t he t a ) )  <  rhom 

L2  =  rho[il]  +  rhom 

LI  =  rho  [  i  1  ]  —  rhom 

12  =  doppler  .integral  (  rho  [  i  1  ]  ,  L2  ,  rhox,  rhoz,  rhom,  theta  [il]) 

II  =  doppler  .integral  (  rho  [  i  1  ]  ,  LI,  rhox,  rhoz,  rhom,  theta  [il]) 
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return_val[il]  =  12  —  II 

L2  =  rho[i2]  +  rhom 

LI  =  — rhox *  cos ( the t a [ i2 ] )  +  rhoz *  sin ( theta [ i2 ] ) 

retur  n  val  [  i2  ]  =  doppler  integral  ( rho  [  i2  ]  ,  L2  ,  rhox  ,  rhoz,  rhom, 
theta  [  i 2  ] ) \ 

—  doppler  integral  ( rho  [  i2  ]  ,  LI,  rhox,  rhoz,  rhom,  theta  [i2  ]) 

L2  =  rhox *  cos ( t he t a  [  i3  ] )  +  rhoz *  sin ( theta  [  i3  ] ) 

LI  =  rho[i3]  —  rhom 

r  etur  n  _val  [  i3  ]  =  doppler  .integral  (rho  [  i3  ]  ,  L2  ,  rhox,  rhoz,  rhom, 
theta  [  i 3  ] ) \ 

—  doppler  .integral  (rho  [  i3  ]  ,  LI,  rhox,  rhoz,  rhom,  theta  [i3  ]) 
return  return.val 

#  classes  for  implementing  functors  to  simplify  some  computations 
class  AntennaPowerPattern(  object  ): 

def  _  _ i  n i  t  _  _  (  se If  ,  gain.db  ,  antenna.depression_angle.deg  , 
beamwidth_deg  ,  patter ntype  =  ’rect’): 

self.gain_db  =  gain_db 
s e 1 f . ante nna .depre s s ion .angle _de g  =  \ 
antenna_depression_angle_deg 
s e If  .  beam width.deg  =  beamwidth.deg 
self,  patter  ntype  =  patter  ntype 

def  _  _c  a  1 1  _  _  (  se  If  ,  phi,  theta): 

if  self. patter  ntype  ==  ’rect’: 
return  sincpattern( 
phi  *180. /pi  , 

theta  *180. /pi  — self,  antenna.depression_angle.deg  , 

self . gain_db  , 

self  .  beamwidth.deg  , 

) 

elif  self. patter ntype  ==  ’radial’: 

theta.a  =  pi / 1 8 0 .* s e If . ante nna  .depre s s ion .angle  _de g 
x  =  cos(theta_a)*cos(theta)*cos(phi)  +  \ 
sin(theta_a)*sin(theta) 
y  =  cos(theta)*sin(phi) 

z  =  —sin ( theta.a )* cos ( theta )* cos ( phi )  +  \ 
cos(theta_a)*sin(theta) 
phi  .prime  =  arctan2(y,  x) 

theta.prime  =  arctan2(z,  sqrt(x**2.  +  y**2.)) 
phi.prime.deg  =  1  8  0 .  /  pi  *  phi  .prime 
theta.prime  _deg  =  1  8  0 .  /  pi  *  t  he  t  a  .prime 
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return  sinc2pattern  ( 
phi  .prime  _deg  , 
theta.prime.deg  , 
self . gain.db  , 
self  .  beamwidth.deg  , 

) 

class  DepressionAngle  (  object  ) : 

def  _  _  i  n  i  t  _  _  (  se  If  ,  altitude): 
self  .  altitude  =  altitude 

def  _  _c  a  1 1  _  _  (  self  ,  tau): 

output  =  zeros  (len(tau)) 

zeroind  =  0. 5  *  SPEED  OF  LIGHT  M  PER  S*  t  au  <  s  e  1  f  .  al  t  i  t  u  d  e 
output[zeroind]  =  pi/2, 
output  [~zeroind]  =  a  r  c  s  i  n  ( 

self  .  altitude  /( 0 . 5  *  SPEED  _OF_LIGHT_M_PER_S*  t  au  ["zeroind  ])) 
return  output 

class  NrcsLookupTable  (  object  ) : 

def  _ init _ (self  ,  f  =  10e9,  polarization=’vv’): 

if  polarization  ==  ’ vv  '  and  8e9  <=  f  <=  12e9: 

data  =  genfromtxt  (  ’  .  /  data  /  sigmaOvv.xband  .  csv  ’  , 
d  e  1  i  m  i  t  e  r  =  ’  ,  ’  ) 

else  : 

raise  Exception/ 

’  invalid  ^polarization^or^frequency^choice’) 
self.x  =  data [ : , 0 ] 
self.y  =  1 0. *  * ( data [ :  , 1 ] / 1 0 . ) 

def  _ _c a  1 1 _ _ ( se If  ,  grazing .angle _rad ) : 

grazing .angle _deg  =  grazing .angle _rad  *  1 80./ pi 
output  =  zeros(len(grazing_angle_deg)) 

badind  =  (  grazing .angle _deg  <  0)  |  (  grazing .angle _deg  >=  90) 

output  [  badind  ]  =  0 

output [~ badind ]  =  interp  ( grazing .angle _deg [~ badind ] , 
self.x,  self.y) 

return  output 

class  FlatEarth DepressionAngle  /  object  ): 

#  also  equal  to  the  grazing  angle 

def  _ init _ (self  ,  altitude): 

self  .  altitude  =  altitude 

def  _  _c  a  1 1  _  _  (  se  If  ,  tau): 

r  =  0.5*SPEED_OF_LIGHT_M_PER_S*tau 
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return.val  =  zeros (tau. shape) 
i  =  r  >=  self,  altitude 

return_val[i]  =  a  rcsin(self.  altitude  /  r  [  i  ] ) 
return_val[~i]  =  pi/2, 
return  return.val 

class  AzimuthAngle  (  object  ) : 

def  _  _  i  n  i  t  _  _  (  se  If  ,  v,  f,  theta): 

if  v  [  1  ]  !=  0:  raise  Exception  (’  v_y  ^must^be  uzero  ’ ) 

self  .  wavelength  =  SPEED _OF_LIGHT_M_PER_S /  f 
self.rho.x  =  2*  v  [  0  ]/ s  e  1  f  .  wavelength 
self.rho.z  =  2*  v  [  2  ]/ s  e  1  f  .  wavelength 
self. theta  =  theta 

def  _ _c a  1 1 _ _ ( self  ,  tau,  rho  ) : 

arccosarg  =  (rho  — se 1 f . rho _z *  sin ( se 1 f . theta ( tau  )))/ ( 
self . rho_x*cos( self . theta (tau  ))) 
badinds  =  abs ( arccosarg )  >  1. 
if  nnz(badinds)  >  0: 

arccosarg [ badinds ]  =  sign ( arccosarg [ badinds ] ) 
return  arccos ( arccosarg  ) 

class  Clutter  Scatter  ingFunction  (object  ): 

def  _  _  i  n  i  t  _  _  (  se  If  ,  f,  sigma_0  ,  v,  phi,  theta,  alpha,  G,  h,  rhom): 
self  .  wavelength  =  SPEED _OF_LIGHT_M_PER_S /  f 
self.sigma_0  =  sigma_0 
if  v [ 1 ]  !=  0: 

raise  Exception  (  ’  v_y  ^must  ^be  ^ zero  ’ ) 
if  v[0]  ==  0  and  v[2]  ==  0: 

raise  Exception(  ’v^must^be^nonzero  ’) 
self.rho.x  =  2*  v  [  0  ]/ s  e  If  .  wavelength 
self.rho_z  =  2*  v  [  2  ]/ s  e  1  f  .  wavelength 
s e If  .  phi  =  phi 
self. theta  =  theta 
self,  alpha  =  alpha 
self  .G  =  G 
s  e  1  f  .  h  =  h 
self.rho.m  =  rho_m 

def  _  _c  a  1 1  _  _  (  self  ,  tau,  rho): 
thetatau  =  self.theta(tau) 
cost h eta  =  cos(thetatau) 
sintheta  =  sin(thetatau) 

il  =  abs(rho  —  s e If . rho _z * s i n t he t a )  <  self . rho_x* costheta  +  \ 
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self. rho_m 

i2  =  0.5  *  SPEED _OF_LIGHT_M_PER_S*  tau  >  self.h 
i  =  i  1  &  i  2 

return.val  =  zeros (tau. shape) 

r  etur  n  _val  [  i  ]  =  s e  If  .  wavelength  **2*SPEED_OF_LIGHT_M_PER_S  / ( 

4 .  *  p  i )  *  *  3 . 

return.val [ i ]  *=  self  .  sigma_0  (  self  .  alpha  ( tau  [  i  ]) )/( 

0.5*  SPEED_OF_LIGHT_M_PER_S  *tau[i])**3. 
r etur n _val [ i ]  *=  s e 1 f ,G(  s e If . phi ( t au [ i ] , rho [ i ] ) , 
thetatau [  i  ] )  *  *  2 . 
return_val[i]  *=  dopple r .slice  ( 

rho  [  i  ]  ,  self.rho.x,  self.rho.z,  self.rho.m,  thetatau [i]) 
return.val [isinf (return.val)]  =  0. 
return_val[isnan(return_val)]  =  0. 

return  return.val 

class  UnmodulatedPulse  Ambiguity  Mag  (object  ) : 
def  _ _ i n i t _ _ ( se If ,  T): 
self .T  =  T 

def  _  _c  a  1 1  _  _  (  se  If  ,  tau,  rho): 
i  =  abs(tau)  <=  self.T 
return.val  =  zeros (tau. shape) 

return_val[i]  =  abs((l  —  abs  (tau[i  ])  /  self.  T)  *  sine  (self.T*  rho  [i]*( 
1— abs ( tau [ i ] ) / self.T))) 
return  return.val 

class  PulseTrainAmbiguityMag(  object  ): 

def  _  _  i  n  i  t  _  _  (  se  If  ,  T,  N,  T_r  ,  chi  .pulse  _mag  ) : 
self.T  =  T 
self  .N  =  N 
self  .  T.r  =  T.r 

self  .  chi.pulse.mag  =  chi.pulse.mag 

def  _  _c  a  1 1  _  _  (  self  ,  tau,  rho): 

return.val  =  zeros (tau. shape) 

i  =  abs  (tau)  <=  self.  N*  self.  T.r 

for  p  in  xrange  ( — (  s  e  1  f  .N— 1) ,  self.N): 

return.val  [  i  ]  +=  1  ./ self  .N*  self  .  chi.pulse.mag  ( 

tau  [i]  —  p*self.T_r  ,  rho[i  ])*  pulsetrainfactor( 
rho  [  i  ]  ,  self.N,  p,  self. T.r) 
return  return.val 


class  PulseTrainPeriodicAmbiguityMag(PulseTrainAmbiguityMag): 
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def  _ _c a  1 1 _ _ ( self  ,  tau  ,  rho  ) : 

tau_shift  =  copy  ( tau )%  self  .  T.r 
i  =  tau  .shift  >  self.Tr/2. 
tau_shift[i]  =  s  e  1  f  .  T  _r  —  tau_shift[i] 

return  1  ./ self  .N*  self  .  chi.pulse.mag  ( tau  shift  ,  rho)  *  \ 
puls  e  .train  .factor  (rho  f  self.N,  p=0,  T_r=  self  .  T_r  ) 

class  PointTargetTimeFreqDistr(  object  ): 

def  _  _  i  n  i  t  _  _  (  se  If  ,  E_x  ,  G,  phi,  theta,  f,  sigma,  tau.bar, 
rho.bar  ,  chi_mag_sq): 

s  e  1  f  .  E_x  =  E_x 
self  .G  =  G 
s e  1  f  .  phi  =  phi 
self. theta  =  theta 

self  .  wavelength  =  SPEED _OF_LIGHT_M_PER_S /  f 

self,  sigma  =  sigma 

self. tau.bar  =  t  a  u  _  b  a  r 

self. rho.bar  =  rho.bar 

s e  1  f  .  chi _mag _sq  =  chi.mag.sq 

def  _  _c  a  1 1  _  _  (  se  If  ,  tau,  rho): 

return  self. E_x *  self. G( self. phi (self. tau.bar  , self. rho.bar), 
self.  theta(self.  tau  .bar))*  *2.  *  \ 
s  e  1  f  .  wavelength  *  *2  .*  s  e  1  f  .  sigma /(( 4  .  *  pi  )**  3 .  *  \ 

(0.5*  SPEED  OF  LIGHT  M  PER  S*  self.  tau_bar)**4.)  *  \ 
self  .  chi.  mag_sq(tau  — self  .  tau.bar  ,  rho  — s  e  If  .  rho.bar) 


#  main  program  code 
def  main  ( ) : 

#  i n puts 
f  =  10e9 
h  =  5e3 

v  =  array ( [300 ,  0 ,  0] ) 
v.wave  =  1.0 
phi.a.deg  =  0. 
theta.a.deg  =  25. 
beamwidth.deg  =  10. 
gain.db  =  20. 

T  =  10e-6 
T.r  =  100e-6 
N  =  8 
E.x  =  N 

fft_pts.min_per.prf  =  256 


Statistical  LTV  System  Model  for  Computing  Clutter  Interference  Power 


27 


r_max  =  50e3 
figsize  =  (5.25,3.5) 

#  Calculated  parameters 

T_s  =  T/10. 

rho _x  =  2*v[0]*  f/ SPEED _OF_LIGHT_M_PER_S 
rho  z  =  2* v [2] *  f  / SPEED  OF  LIGHT  M  PER  S 
rho.wave  =  2*v_wave*f  /  SPEED  _OF_LIGHT_M_PER_S 
print  ’ rho.wave ,  rho_wave 

rho_h  =  rho _x  *  cosd  ( t he  t  a _a _de g  ) *  cosd  (  phi  _a_de g  +  beamwidth_deg /2  . )  +  \ 
rho_z*sind(theta_a_deg  ) 

deltarho  =  rho_x*cosd(theta_a_deg)*sind(phi_a_deg)*( 
beamwidth.deg  /2.*pi/180.) 

print  ’  differential  ^doppler  ^  be  am  width  ^(onesided)  ’  ,  deltarho 

r_a  =  h*cscd(theta_a_deg) 

tau_a  =  2*r_a  /SPEED  OF  LIGHT  M  PER  S 

rho_a  =  rho_x*cosd(theta_a_deg)*cosd(phi_a_deg)  +  rho_z*sind( 
theta_a_deg  ) 

print  ’tau_a^=^’,  tau_a*le6 
print  ’rho_a^=^’,  rho_a 

print  ’  doppler  ^beamwidth  ^  (  onesided  )  ’  ,  abs(rho_h  —  rho_a) 

if  abs(rho_h  —  rho_a)  <  rho_wave : 

raise  Exception!  ’  antenna^ approximation^not^valid  ’) 
tau.max  =  2*r_max  / SPEED _OF_LIGHT_M_PER_S 
if  tau_a  >  tau_max: 

raise  Exception!  ’  target^is  ^outside  ^max^  simulation^  range  ’) 

B  =  rho  h  +  rho_wave  +  1/T_r  #  clutter  bandwidth 

#  functors  for  doing  computations 
G  =  AntennaPowerPattern  ( 

gain_db  =  gain_db  , 

antenna_depression_angle_deg=theta_a_deg  , 
beamwidth_deg  =  beamwidth_deg  , 
patter ntype  =  ’  rect  ’ 

) 

theta  =  FlatEarthDepressionAngle  ( 
altitude  =h, 

) 

alpha  =  theta 

phi  =  AzimuthAngle  ( v ,  f,  theta) 
sigma_0  =  NrcsLookupTable  (  f ) 

def  G_rot(x,y):  return  G(x— abs  (  phi  _a_deg  *  pi  /  1  8  0 . )  ,  y) 

C_eta  =  ClutterScatteringFunction  (f ,  sigma.O  ,  v,  phi,  theta, 
alpha,  G_rot  ,  h,  rho_wave) 


#  time  and  frequency  axis  computations 
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if  N<  fft_pts_min_per_prf: 

fft_pad_factor  =  fft_pts_min_per_prf/N 
else  : 

fft_pad_factor  =  1 
drho  =  l./(N*T_r)/fft_pad_factor 
print  ’ drho : ’ ,  drho 
print  ’  PRI  ^=J7cg  ^us  ’  %(T_r  *  1  e6 ) 
print  ’PRF.=  J7cg  JcHz’%(le-3/T_r) 
if  drho  >  rho_wave  :  raise  Exception  ( 

’  frequency  ^  re  s olution  ds  doo  dow  ’ ) 
print  ’  tau.max  :  ’  ,  tau_max*le6,  ’us’ 

Npri  =  int  (  c  e  i  1  ( tau.max  /  T_r  ) ) 
if  (Npri%2)  ==  0:  Npri  +=  1 

#  pulse  train  ambiguity  computation 
chi_pulse_mag  =  UnmodulatedPulseAmbiguityMag  (T) 
print  ’  computing  ^pulsedrain^ambiguity  ’ 

chi  train  mag  =  PulseTrainPeriodic  AmbiguityMag  (T ,  N,  T  r, 
chi_pulse_mag  ) 

def  chi  -train _mag _s q  (  ttt  ,  rrr):  return  chi_train_mag  (  ttt  ,  rrr)**2. 

print  ’  done  ^computing  ^pulse  strain  ^ambiguity  ’ 

Ntau  =  int  ( round  (  T_r  /  T_s  ) ) 
tau_chi_lin  =  arange  ( Ntau )  *  T_s 
rho_chi_lin  =  arange(— B,  B,  drho) 

Nrho  =  len(rhochilin) 

tau_chi  ,  rho_chi  =  meshgrid  (  tau  _chi  _lin  ,  rho_chi_lin) 

#  plotting  pulse  train  ambiguity 
print  ’  plotting  ^pulse^  tr  ain^ambiguity  ’ 
tt  =  ( tau_chi  —  T_r  12  . ) 

rr  =  (  rho_chi  —  0.5/T_r ) 

xx  =  db(chi_train_mag_sq(tt  ,  r r  ) ) 

plot -image  (  1 1  *  1  e6  ,  rr/le3,  xx  —  xx.max(),  figsize) 
xlabel  (r  '  $\ tau$  ^ ( $\mu$s )  ’ ) 
y lab el ( r ’ $\rho$  ^(kHz) ’ ) 

title  (r’$|\chi_  {NT  }(\tau  ,\rho)|"2$u(dB^normalizeduto^peak)  ’) 

ylim(—  0.5/T_r/le3  ,  0.5/T_r/le3) 

dim  (  —20  ,0) 

tight-layout  () 

print  ’ done ^ p 1 o 1 1  i  ng ’ 

savefig  (  ’  fig.trainambig  .pdf’) 

#  repmatting  ambi guity  function 
print  ’  pri  u  ratio  :’  ,  Npri 

c  hi -train  _mag  _sq  _repmat  =  tile  (  chi -train  _mag  _sq  ( tau.chi  ,  rho.chi), 
(1,  Npri)) 
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tau  =  zeros  (  c  hi  .train  _mag  _sq  .repmat  .  shape  ) 

rho  =  zeros  (  c  hi  .train  _mag  _sq  .repmat  .  shape  ) 

for  i  in  xrange  (  Npri  ) : 
for  j  in  xrange ( 1 ) : 

tau  [  j  *Nrho  :  ( j  +l)*Nrho  ,  i  *Ntau  :  (  i  + 1  )*  Ntau  ]  =  tau_chi  +  i*T_r 

rho  [  j  *Nrho  :  ( j  +l)*Nrho  ,  i  *Ntau  :  (  i  + 1  )*  Ntau  ]  =  rho_chi  +  j/T_r 

rho  — =  (l//2)/T_r  #  centering  Doppler  axis  on  0  Hz 

#  C-eta  computation 
print  ’  C _et a  ^ c  ale  ’ 

i  =  abs  ( rho  _chi  _lin  )  <=  rho_h  +  rho_wave  +  0.5/T_r 
C_eta_mat  =  C_eta(tau[i  ,:]  ,rho[i  ,:]) 
print  ’  done  u  with  _eta  ’ 

#  output  power  distribution  computation 
print  ’  convolving  ^ ambiguity  ’ 
cetarhomax  =  max(  r  ho  _c  hi  _1  in  [  i  ] ) 
cetarhomin  =  min(  rho  chi  lin  [  i  ]) 
chirhomax  =  max(  r h o  _c  hi  _1  i  n  ) 
chirhomin  =  min (  r h o  _c  hi  _1  i  n  ) 

rhomax  =  cetarhomax  +  chirhomax 
r  h  o  m  i  n  =  cetarhomin  +  chirhomin 
minind  =  rhomin/drho 
maxind  =  rhomax  /drho 

rhonew  =  arange  ( minind  ,  maxind  +  l)*drho 
i  =  find((0  <=  rhonew)  &  (rhonew  <  l./T_r)) 
tau  P  =  arange  (  Ntau  )*  T_s 

rho_P  =  arange  (  f  f  t  _p  t  s  _mi  n  _p  e r  _prf  )*  drho 
tau_P  ,  rho_P  =  meshgrid  ( tau_P  ,  rho_P) 

P  =  E_x*  circconv2 _tauaxis  (  C_eta_mat  ,  chi_train_mag_sq_repmat 
) .  real  *T_s*drho 
P  =  P  [  i  , :  Ntau] 

print  ’  done  ^convolving  ^ambiguity  ’ 

plot_image  ( tau_P  *le6  ,  rho_P/le3,  db  (P/P .  max  ( ) )  ,  figsize) 
xlabel(r’$\tau$u  ( $\mu$s )  ’ ) 
ylabel(r'$\rho$  ^(kHz) ’  ) 

title  (  r  ’  $P_c  (\  tau  ,\rho)$^(dB^normalized^to^peak)  ’) 
clim(—  20,  0) 
tight_layout  () 

savefig  (  ’  fig.outputpdmap  .  pdf  ’  ) 

#  point  target  return 
sigma  =  10000 

P.prime  =  PointTarge  tTimeFreqDis  tr  ( E_x  ,  G_rot  ,  phi,  theta,  f, 
sigma,  tau_a,  rho_a  ,  c  hi  _tr  ai  n  _mag  _sq  ) 
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Pp  =  P .prime  ( tau_P  ,  rho_P) 

plot.image  ( tau_P  *le6  ,  rho_P/le3,  db  ( Pp /Pp  .  max  () )  ,  figsize) 
xlabel(r’$\tau$^  ( $\mu$s )  ’ ) 
ylabel(r’ $\rho$  ^(kHz) ’ ) 

title  (r”$P_t(\tau  ,  \  rho  )$  u  (dB  ^  normalized  „  to  ^peak )”  ) 
clim(— 20,  0) 
tight.layout  () 

savefig  (  ’  fig.pointtarget  .  pdf’  ) 

#  calculating  signal— to  — clutter  ratio 
tau _a_aliased  =  tau_a 
rho_a_aliased  =  rho_a 

while  t  a  u  _a  _al  i  as  e  d  >  T_r:  tau  _a_aliased  — =  T_r 
while  r  ho  _a  a  1  i  a  se  d  >  l./T  r:  rho  a  aliased  — =  l./T  r 
i  =  argmin  (  abs  ( tau  _P  [0  , —  t  a u  _a  _al i  a s e  d  ) ) 
j  =  argmin  ( abs  ( rho  P  [:  ,0]  —  rho  a  alias  ed  ) ) 
print  ’ tau _a_aliased „= ’  ,  t a u _a _al i a s e d 
print  ’rhoaaliased^=’,  rhoaaliased 

sigpwr  =  P.prime  (  array  ([  tau  _a  _  a  1  i  as  ed  ])  ,  array  ([  rho  _a_ali  ased  ] )) 

clutter  pwr  =  P[j,i] 

print  ’  db  (  sigpwr  )  ,  db(sigpwr) 

print  ’db(clutterpwr)^=’  ,  db  (clutter  pwr) 

print  ’  SCR.discrete  ,  db  (  sigpwr  /  clutterpwr  )  ,  ’dB’ 

#  compare  Richards  ’’Fundamentals  of  Radar  Signal  Processing”  p.48 
if  TV  t  au  _a  *  tand  ( t  he  t  a  _a  _deg  )  >  beamwidth  deg  *  pi  /  1  8  0 . : 

print  ’beanie limited  ^according  ^to  ^Richards  ’ 
else  : 

print  'pulse  ^limited  ^according  ^to  ^Richards  ’ 

R  =  0.5*SPEED_OF_LIGHTJVLPER_S*tau_a 
DeltaR  =  0 .5  *  SPEED _OF_LIGHT_M_PER_S*T 
SCR_richards_blimited  =  sigma*  sind  ( theta_a_deg  )/( 

R*  *2 .  * (  beamwidth.deg  *  pi  /180.)**2.*  sigma.O  ( 
array  ([  theta_a_deg*pi  /  I  80.]))) 

SCR_richards_plimited  =  sigma*cosd  ( theta_a_deg  )/( 

R*  sigma.O  (  arr  ay  ( [  the  ta  _a  _deg  *  pi  / 1  8  0 .  ] ) )  *  \ 

DeltaR*  be  am  width  _deg*  pi  /  1  80.) 

print  ’  SCR_richards_blimited  ,  db(SCR_richards_blimited),  ’dB’ 
print  ’  SCR  richards  plimited  ,  db(  SCR  richards  plimited  )  ,  ’dB’ 

#  plotting  SCR 
SCR  =  Pp/P 

plot  image  ( tau  P  *  1  e6  ,  rho  P/le3,  db(SCR),  figsize) 
xlabel(r’$\tau$^  ( $\mu$s )  ’ ) 
ylabel(r'$\rho$  ..(kHz) ’ ) 
title  (r”$P_t/P_c$^(dB)”) 
c  1  i  m  (  — 10,  10) 
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t  ig  ht  _1  ay  o  u  t  () 
savefig  (  ’  fig  _scr  .pdf’) 


if  name  ==  ’  main  ’  : 
main  () 


